Exchange-Correlation Hole of a Generalized Gradient Approximation 

for Solids and Surfaces 



Lucian A. Constantin 1 , John P. Perdew 1 , and J. M. Pitarke 2,3 
1 Department of Physics and Quantum Theory Group, Tulane University, New Orleans, LA 70118 
2 CIC nanoGUNE Consolider, Mikeletegi Pasealekua 56, E-20009 Donostia, Basque Country 
3 Materia Kondentsatuaren Fisika Saila (UPV/EHU), 
DIPC, and Centro Fisica Materiales CSIC-UPV/EHU, 
644 Posta kutxatila, E-48080 Bilbo, Basque Country 
(Dated: January 25, 2009) 

We propose a generalized gradient approximation (GGA) for the angle- and system-averaged 
exchange-correlation hole of a many-electron system. This hole, which satisfies known exact con- 
straints, recovers the PBEsol (Perdew-Burke-Ernzerhof for solids) exchange-correlation energy func- 
tional, a GGA that accurately describes the equilibrium properties of densely packed solids and their 
surfaces. We find that our PBEsol exchange-correlation hole describes the wavevector analysis of 
the jellium exchange-correlation surface energy in agreement with a sophisticated time-dependent 
density-functional calculation (whose three-dimensional wavevector analysis we report here). 
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I. INTRODUCTION 

In the Kohn-Sham (KS) density functional theory^ for 
the ground-state energy of a many-electron system, only 
the exchange-correlation (xc) energy £J xc [n] has to be ap- 
proximated. The exact xc energy of an arbitrary inho- 
mogeneous system of density n(r), which incorporates all 
the quantum many-body effects beyond the Hartree ap- 
proximation, can be obtained from the spherical average 
n xc (r, u) of the coupling-constant averaged xc hole den- 
sity n xc (r, r') at r' around an electron at r as follows^ 3 - 
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drn(r) s xc [n] (r), 
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where e.;; C [n](r) is the xc energy per particle at point r 

1 f°° 1 
£ X c[n](r) = - / du 4:Tru 2 -n xc (r,u) 
£ Jo u 

= 4 Jo°° dk io°° du « 2 H^T Mr, u), 

with 

n xc (r,u) = -!- J dSln xc (r,r'), 
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dfi being a differential solid angle around the direction 
of u = r' - r, and k representing the magnitude of the 
wavevector. (Unless otherwise stated, atomic units are 
used throughout, i.e., e 2 = % — m e = 1.) 

The "Jacob's ladder" classification of the widely- 
used ground-state density-functional approximations for 
£Ja; C [n] and n xc (r, u) has three complete non-empirical 
rungs: the local-spin-density approximation (LSD A), 1 
the generalized-gradient approximation (GGA),— ^ and 
the meta-GGA»£i£ Due to its simplicity and accuracy, one 
of the most commonly used xc density-functional approx- 
imation in solid-state physics and quantum chemistry cal- 
culations is nowadays the semilocal PBE-GGA4 



Recent work-2 has shown, however, that the exchange 
density-functional approximations should recover, in the 
limit of slowly-varying densities, the universal second- 
order gradient-expansion (GE) approximation of the ex- 
change energy, 
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where £ x m f is the exchange energy per particle of the 
uniform electron gas, /i x E = 10/81 is the GE exchange 
coefficient,— and s — \'Vn\/(2kpn) is the reduced den- 
sity gradient which measures the variation of the elec- 
tron density over a Fermi wavelength Xp = 2tt /kp, with 
kp = (Sir^ri) 1 ^ representing the magnitude of the lo- 
cal Fermi wavevector. Recovery of the correct second- 
order gradient expansion for correlation^ in the slowly- 
varying limit is much less important for the construction 
of density-functional approximations. (See Table 1 of 
Ref.[l). " 

A GGA, which has as ingredients only the spin densi- 
ties ri| and 7i j, and their gradients Vrij and Vn^, cannot 
recover, in the limit of slowly varying densities, the GE 
approximation of the exchange energy and at the same 
time be accurate for atoms^^ The semilocal PBE has 
the correct correlation GE coefficient in the high-density 
limit, and is accurate for atoms, but its exchange GE 
coefficient is almost twice as large as the exact coeffi- 
cient, i.e., [i BBE w 2p^ E . Because of this^ PBE over- 
estimates the equilibrium lattice constants of solids and 
yields surface energies that are too low. 

Following the ideas of Ref. ©, PBEsol (PBE for solids) 
was constructed^. PBEsol is a GGA that has the same 
form as PBE but restores the density-gradient expan- 
sion for exchange by replacing /j, x be = 0.2195 with 
fi PBEso1 — n x E ■ By fitting the jellium xc surface en- 
ergies (as had been done previously in Ref. 13J in the 
construction of a GGA which relies on the Airy-gas 
approximation^), the PBEsol correlation GE coefficient 



2 



was set to (i PBEso1 = 0.046. (For PBE, fi PBE = 0.0667). 
Thus, PBEsol can easily be applied in solid-state calcu- 
lations (just by changing the coefficients in a PBE code) 
and yields good equilibrium lattice constants and jellium 
surface energies^ Several other applications of PBEsol 
have already proved the accuracy of this GGA for solids. 
In particular, PBEsol considerably improves the struc- 
ture of gold clusters^ and works better than PBE for 
isomerization energies and isodesmic stabilization ener- 
gies of hydrocarbon molecules^ PBEsol also describes 
ferro- and anti-ferro-electric ABO3 crystals^ much bet- 
ter than LSDA or PBE-GGA. 

In this paper, we first construct the PBEsol angle- 
averaged xc hole density h PBEsol (r,u). A nonempirical 
derivation of the PBE xc hole was reported in Ref. |5j, 
starting from the second-order density-gradient expan- 
sion of the xc hole and cutting off the spurious large- 
it contributions to satisfy exact constraints according to 
which (i) the exchange-hole density must be negative, (ii) 
the exchange hole must integrate to -1, and (iii) the cor- 
relation hole must integrate to zero. Later on, a fully 
smoothed analytic model was constructed for the PBE 
exchange hole^ Our construction of the PBEsol xc hole 
begins with and appropriately modifies the sharp cutoff 
correlation hole of Ref. [E] and the smooth exchange hole 
of Ref. @. It should be recalled that, because of an inte- 
gration by parts that occurs in the underlying gradient 
expansion, a GGA hole is meaningful only after averag- 
ing over the electron density n(r) (as in our tests and 
applications), and this system- averaging itself smooths 
sharp cutoffs. 

Finally, we use our PBEsol xc hole to carry out a 
three-dimensional (3D) wavevector analysis of the jellium 
xc surface energy. This wavevector analysis was carried 
out in Ref. [H in the random-phase approximation (RPA). 
Here, we go beyond the RPA in the framework of time- 
dependent density- functional theory (TDDFT), and we 
compare these calculations with the results we obtain 
from our PBEsol xc hole density. 

The paper is organized as follows. In Sec. II, 
we present the PBEsol angle-averaged xc-hole density 
Sj C B£so! (r,w). In section III, we perform the wavevec- 
tor analysis of the jellium xc surface energy In Sec. IV, 
we summarize our conclusions. 



II. PBESOL-GGA ANGLE-AVERAGED 
EXCHANGE-CORRELATION HOLE 

We assume here that the PBEsol correlation energy 
can be constructed from a gradient expansion for the 
correlation hole in much the same way that the PBE 
correlation energy was so constructed^. The GGA angle- 
averaged correlation hole is£ 

n GGA (r Sl C, t, v) = 5 fc 2 L4 c (r s , (, v)+t 2 B c (r s , £ v)}6(v c -v), 

(5) 

where r s — (9tt /A) 1 / 3 /Isf is a local density parameter, 
£ = (n-f — ^x)/( n T + n i) IS the relative spin polariza- 



tion, (f) = [(1 + C) 2/3 + (1 - C) 2/3 ]/2 is a spin-scaling 
factor, v — 4>k s u with k s = {Akp/i:) 1 / 2 is the reduced 
electron-electron separation on the scale of the screen- 
ing length, and t — | Vn|/(20fc s n) is the reduced density 
gradient measuring the variation of the electron density 
over the screening length. The sharp cutoff v c is found 
such that Eq. (O satisfies the correlation hole sum rule 
/ dr ra c (r,r') = 0. (j> 5 k 2 A c (r s , (, v) is the LSDA correla- 
tion hole given by Eq. (45) of Ref. 5, and the gradient 
correction to the correlation hole is given by the following 



expression:" 



B c (r s ,C,v) = B^(v)[l- 



+ /3(r s ,C)uV 



(6) 



where B L (v) is the RPA nonoscillating long-range 

contribution given by Eq. (49) of Ref. [3 p(r s ,Q — 

7t/cf(0.305 — 0.136£ 2 )/4</> 4 measures where the short- 
range contribution vanishes, and 
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is constructed so that the second-order gradient expan- 
sion of the PBEsol correlation energy is recovered. Here 
Ei(y) = U eV dt e~* ft is between and 1, and p. GGA is 
the GGA gradient coefficient in the slowly-varying limit 
(fi PBE for PBE and n PBEso1 for PBE-sol). 

In Fig. [TJ we show PBE and PBEsol versions of 
B c (r s ,(,v) versus v, for r s = 2. Both PBE and PBEsol 
recover the correct RPA-like behavior [B^ («)] at large 
v, and they both show the same £ behavior; because 
^PBEsol < n PBE : however, at intermediate values of v 
the PBEsol gradient correction to the correlation hole is 
substantially smaller than the PBE one. The gradient 
correction B c (r s ,(, v) of Eq. ([6]) can be negative at small 
v and small r s ; however, because of the energy sum rule 
both f^°du uB c (r s ,£,v) and du u 2 B c (r s ,(,v) are 
positive, which ensures that the cutoff procedure is cor- 
rect and for every value of r s , £, and t there is a v c such 
that Eq. ([5]) satisfies the correlation-hole sum rule. 

The exchange energy and exchange-hole density for 
a spin-polarized system may be evaluated from their 
spin-unpolarizcd counterparts by using the spin-scaling 
relationship 



E x [n r ,n l } = -{E x [2n^+E x [2 ni \} (8) 



and 



»)=£ 



n{r) 



n x [2n a }{Y, r + u); (9) 



thus, we only need to consider the spin-unpolarized sys- 
tem. As in the case of the analytical PBE exchange hole 
of Ref. @, we choose the following ansatz for the nonoscil- 
latory dimensionless exchange-hole shape: 
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FIG. 1: Gradient correction to the correlation hole, 
B c (r a ,C,v), versus v for r s = 2. PBE and PBEsol are com- 
pared here. The solid line represents the Langreth-Mehl (LM) 
RPA contribution, which should be recovered at large v. 




FIG. 2: The exponent H{s) of Eq. (fT0|) versus the reduced 
gradient s. The solid line represents the numerical solution 
of the implicit equation for TL{s) [Eq. (A4) of Ref. [6j. The 
dashed line represents the fit of Eq. (|14|l . 



where s is the reduced density gradient for exchange. 
When s = 0, Eq. (TOD recovers^ J LSDA for A = 
1.0161144, B = -0.37170836, C = -0.077215461, V = 
0.57786348, and £ = -0.051955731. The functions F(s), 
Q(s), and H(s) are found in such a way that the energy 
and exchange-hole sum rules arc satisfied: 
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dy yj 



PBEsol , 



s,y) 



-F 



PBEsol 



(s) 



and 



f / d V V>J PB *«*{a tV ) = -l, 



(11) 



(12) 



and also the small-it behavior of the exchange hole is 
recovered by£: 



T{s) = 6A75H{s) + 0.4797. 



(13) 



Here F PBEsol (s) is the PBEsol enhancement factor^ 
The integrals of Eqs. (1 1 1 [) and (fl"2f can be solved 
analytically^ and Eqs. (fl~T]) - ([i"3"j) reduce (by substitution) 
to an implicit equation for H [Eq. ( A4) of Ref. @] . We 
have solved this equation for PBEsol; the numerical so- 
lution that we have found for H(s) (see Fig. [2]) can be 
fitted to the following analytic expression: 



H(s) 



a 2 s 



1 + a 3 s 4 + 04s 6 



(14) 



where ai = 0.00018855, a 2 = 0.00741358, a 3 = 
0.05687256, and a 4 = 0.00675093. For s > 8.5, as oc- 
curs in the tail of an atom or molecule where the electron 
density is negligible, the implicit equation for Ti.(s) does 
not have a solution (as in the PBE case^) so we reset s to 
s = 8.5. Recently, Henderson et al? 1 constructed a GGA 
exchange hole that eliminates this unphysical large-s be- 
havior by using some ideas from the meta-GGA holei 
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FIG. 3: Dimensionless exchange hole shape J PBEsoL (s, y) [see 
Eq. (|10[) ] versus y = tu, for s between and 3 in steps of 
0.5. For comparison with J FBE (s,y), see Fig. 2 of Ref. S 



When s = 0, J* 



l {s,y) yields J 



In Fig. [3j we plot the dimensionless exchange hole 
shape J PBEsol {s,y) [using the analytical fit of Eq. (fHj) ] 
versus y = kpu for several values of the reduced gradi- 
ent s. Our J PBEsol {s,y) looks similar to the J PBE (s,y) 
of Ref. 0, but J PBE (s,y) is deeper because [i PBE — 
0.2195 > fi PBEso1 = 0.1235. 

Finally, we look at the xc enhancement factor, which 
displays the nonlocality: 22 



F 



GGA 



^GGA 
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(n) 



(15) 



e umf ^ {j em g exchange energy per particle of a spin- 
unpolarized uniform electron gas. For a spin-unpolarized 
system in the high-density limit (r s — > 0) the exchange 
energy is dominant and Eq. (15) defines the exchange 
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FIG. 4: The PBEsol enhancement factor F X c for the spin- 
unpolarized case (f = 0), as a function of the reduced gra- 
dient s for several values of r 3 . The lines represent the en- 
hancement factor obtained from the PBEsol xc energy func- 
tional of Ref. Q3- The dots represent the enhancement fac- 
tor obtained from our PBEsol angle-averaged xc hole density 
through Eq. ©. 



FIG. 5: The PBEsol enhancement factor F xc for the fully- 
spin-polarized case (£ = 1), as a function of the reduced gra- 
dient s for several values of r 3 . The lines represent the en- 
hancement factor obtained from the PBEsol xc energy func- 
tional of Ref. Q3- The dots represent the enhancement fac- 
tor obtained from our PBEsol angle-averaged xc hole density 
through Eq. ©. 



enhancement factor F X JGA = e GGA (n, Wn)/e x mf (n). In 
Figs. [4] and [5j we show the PBEsol enhancement factor 
for a spin-unpolarized system, F x D c BEsol (r s , £ = 0, s), and 
for a fully-spin-polarized system, F x BEsol (r s ,Q = 1, s), 
versus s for several values of r s , F x BEso1 is calculated ei- 
ther (i) from the analytic expression of e x BEso1 reported 
in Ref. [12] or (ii) from our PBEsol angle-averaged xc hole 
density through Eq. @. Overall, these calculations of 
F x BEso1 agree well with each other, confirming the as- 
sumption made at the beginning of this section; only for 
r s > 10 (when the electron density is very small) and 
s w 1.5 is the error introduced by the second procedure 
significant—. The analytic fit for Ti.(s) used to construct 
our PBEsol exchange hole does not exactly reproduce the 
PBEsol enhancement factor, but the difference is small 
as shown in Fig. [JJ At this point, we also note that 
the parametrization^tS. of H(r s ,C,,t) entering the ana- 
lytic expression of e PBEso1 reported in Ref. [l2| does not 
reproduce exactly the real-space cutoff results, as shown 
in Figs. 7 and 8 of Ref. H. 



III. WAVEVECTOR ANALYSIS OF THE 
JELLIUM XC SURFACE ENERGY 

The xc surface energy, <j xc , can be defined as the xc 
energy cost per per unit area to create a planar surface 
by cutting the bulk. In a jellium model, in which the 
electron system is translationally invariant in the plane 
of the surface, and assuming the surface to be normal to 
the z-axis, the surface energy can be written as follows 3 



where 



(Jl\ 



(16) 



J, Z' + OO 

j xc (k) = 2 — / dzn(z)b xc (k,z) (17) 

7T J_rr, 



is the wavevector-resolved xc surface energy, and 



sinfcit 
ku 



n xc (z,u)-n^(u)] 



b xc (k,z) = 4-7T / duu' 
Jo 

(18) 

Equations ([TB])-([T5]) comprise an angle-averaged three- 
dimensional wavevector analysis of the surface xc energy 
into contributions from density fluctuations of various 
wavevectors k, following from the Fourier transform of 
the Coulomb interaction in Eq. 12]). In these and sub- 
sequent equations, kp = (iit^n) 1 / 3, is the bulk (not the 
local) Fermi wavevector, and r s is the bulk (not the local) 
density parameter. 

The exact low-wavevector limit of a xc is known to be^ 



kp 1 
7xc(k -> 0) = — (u s - -u p )k, 

iTT Z 



(19) 



where u p = ('inn) 1 / 2 and u s = u p /V2 are the bulk- 
and surface-plasmon energies, and n is the bulk density. 
Eq. (|19[) was used in the wavevector-interpolation ap- 
proach reported in Refs.0,0, and[2f| and it was naturally 
recovered by the RPA approach reported in Ref. H. 

Taking into account that k = (k||,k z ), k|| being a 
wavevector parallel to the surface, Eq. (fl"8|) can be ex- 
pressed as follows^ 



b xc (k,z) 



dkz 
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fa* e ik*(*-z') 
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(20) 



In the case of RPA and TDDFT calculations, we 
use the fluctuation-dissipation theorem 2 ^ 26 ' 27 to de- 
rive n xc (k\\;z,z') and n"™*^(fc) from the coupling- 
constant dependent density-response functions ^ (t) 
and XA (fell ; z, z'), as follows 3 -^ 
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x x\( z , z> ; k\\, iu) - n(z)8(z - z')] (22) 



Xa™^ (k, u) and X\( z i z '> k\\ j being 3D and 2D Fourier 
transforms of the corresponding density-response func- 
tion xa(i", r'; w). In the framework of TDDFT (our 
benchmark for this work), the density- response function 
XA( r , r '; w ) satisfies a Dyson-like equation of the form 2 ^ 



XA(r,r';w) = xo(r,r';cj) + / dri dr 2 Xo(r,i*i; 
A 



ri - r 2 



+ /xc,AH( r b r 2^) XA(r2,r';w), (23) 



where xo( r J r ';^) is the density-response function of 
non-interacting KS electrons (which is exactly express- 
ible in terms of KS orbitals^) and f xc ,x[n](r, r'; w) is 
the unknown A-dependent dynamic xc kernel. When 
/ ICj ^[n](r,r';w) is taken to be zero, Eq. (|2"3"]) reduces to 
the RPA density-response function. If the interacting 
density response function XA(r,r';w) is replaced by the 
noninteracting KS density-response function xo(r, r';u;), 
then Eqs. (f2"Tj) and (|2"2"|) yield their exchange-only coun- 
terparts. 

In the calculations presented below, we have consid- 
ered, as in Ref. ,3, a jellium slab of background thickness 
a = 2.23Ai? (where A^ = 2ir/kp) and bulk parameter 
r s = 2.07. This slab corresponds to about four atomic 
layers of Al(100). 

For the GGA calculations of J xc (k), the function 
b xc (k, z) entering Eq. (TP?) is taken from Eq. (fl%|) with 
the xc-hole densities calculated as reported in (i) Ref. [2CJ 
for n™ l/ (u), (ii) Refs.ianddfor n^ E (z,u), and (iii) 
Section II above for n™ Esol {z, u). 

For the exact- exchange, exact-RPA, and TDDFT 
calculations of J xc (k), the function b xc (k,z) entering 
Eq. (fT7"|) is taken from Eq. (f20|) with the xc hole den- 
sities calculated from Eqs. (f2"Tj) and (|22p . In the case of 
the TDDFT calculations, we use the accurate static xc 
kernel reported and used in Ref. [3l|. This kernel, which 
is based on a parametrization 32 of the diffusion Monte 
Carlo (DMC) calculations reported in Ref. [H for the uni- 
form electron gas, was constructed for jellium surfaces 
where neglect of the w-dependence does not introduce 
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FIG. 6: PBE, PBEsol, and exact wavevector-resolved ex- 
change surface energies 7 x (fe), versus k/2kp, for a jellium 
slab of thickness a = 2.23Af and r s = 2.07. The semilocal 
PBE and PBEsol calculations have been performed from non- 
oscillatory parametrizations of the dimensionless exchange- 
hole shapes J PBE and j PBEsoL f respectively. The area under 
each curve represents the corresponding exchange surface en- 
ergy: <rf SB = 2164 erg/cm 2 , crf SBsoi = 2424 erg/cm 2 , and 



~ exact - 2348 erg/cm 2 



significant errors, and is expected to yield exact results 
in the limits of small and large wavevectors. Our numer- 
ical scheme was described in detail in Ref. 0, where only 
RPA calculations were reported. 

The wavevector-resolved exact-exchange surface en- 
ergy 7a, (k) is shown in Fig. [5] Figures[7][5]and Fig.[5]show, 
respectively, the wavevector-resolved xc and correlation- 
only surface energies j xc (k) and J c (k)- Figure [5] shows 
that J x BEsol (k) improves over PBE, as expected, and 
is close to the exact J x {k) for intermediate values of the 
wavevector. At large values of the wavevector, both PBE 
and PBEsol correctly recover the nonoscillatory LSDA 
(see Figs. 1 and 2 of Ref. 0); differences between this 
nonoscillatory PBEsol (and also LSDA and PBE) and 
the exact j x (k) at these large values of k are due to inac- 
curacy of the nonoscillatory model employed in Eq. (jTUJ) 
near k = 2kp. 

In Fig. [71 we compare the wavevector-resolved exact- 
RPA surface energy (as reported in Ref. d) with its 
TDDFT counterpart (which we have not reported else- 
where, and which required three months of computa- 
tion). In the long- wavelength limit (k — > 0), both 
RPA and TDDFT calculations approach the exact low- 
wavevector limit of Eq. (|19[) . In the large- A: limit, 
where RPA is wrong, our TDDFT approach (which re- 
produces accurately the xc energy of the uniform elec- 
tron gas) is expected to be accurate. Furthermore, the 
uniform-gas-based isotropic xc kernel that we use in our 
TDDFT calculation has been shown recently to yield 
essentially the same two-dimensional (2D) wavevector 
analysis as a more sophisticated high-level correlated ap- 
proach (the inhomogeneous Singwi-Tosi-Land-Sjolander 
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FIG. 7: Exact-RPA and benchmark TDDFT wavevector- 
resolved xc surface energies ^ xc (k), versus k/21tF, for a jellium 
slab of thickness a — 2.23Af and r s = 2.07. The area under 
each curve represents the corresponding xc surface energy: 
o PPA = 3091 erg/cm 2 and a PEDFT = 3090 erg/cm 2 . The 
straight dotted line represents the universal low-wavevector 
limit of Eq. iQSJ. 



FIG. 9: PBE, PBEsol, exact-RPA, and benchmark TDDFT 
wavevector-resolved correlation surface energies r y c (k), ver- 
sus k/2kF, for a jellium slab of thickness a = 2.23Af and 
r s = 2.07. The area under each curve represents the corre- 
sponding correlation surface energy: a PBE = 720 erg/cm 2 , 
a PBEsoi = 604 erg/cm 2 , a PPA = 743 erg/cm 2 , and 

a TDDFT = 742 erg/cm 2. 
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FIG. 8: PBE, PBEsol, and benchmark TDDFT wavevector- 
resolved xc surface energies 'yxcik), versus k/2kF, for a jellium 
slab of thickness a = 2.23Af and r s — 2.07. The semilocal 
PBE and PBEsol calculations have been performed from non- 
oscillatory parametrizations of the dimensionless exchange- 
hole shapes J PBE and j PBEso1 ^ respectively. The area un- 
der each curve represents the corresponding xc surface en- 
ergy: a PBE = 2885 erg/cm 2 , a P c BEsoL = 3027 erg/cm 2 , and 
a FEDFT — 3090 erg/cm 2 . The straight dotted line represents 
the universal low-wavevector limit of Eq. (|19p . 



method) which does not use an isotropic kernel derived 
from the uniform gasi^i Hence, we take the TDDFT 
wavevector-resolved surface energy represented in Fig. 
by a solid line as the benchmark curve against which we 
compare various GGA's. 

Figure [5] shows our wavevector-resolved TDDFT sur- 
face energy together with its PBE and PBEsol counter- 



where it matches the exact initial slope of Eq. (jTt 
At intermediate wavevectors, "f xl ? Es ° l IS ver y close to 
our benchmark TDDFT calculation. As in the case 
of wavevector-resolved exchange-only surface energies, 
PBE and PBEsol calculations correctly recover the 
nonoscillatory LSDA and differ from the more accu- 
rate TDDFT calculation due to the inaccuracy of the 
nonoscillatory model of the exchange-hole shape [see, 
e.g., Eq. HUD]. Figure [5] shows that j xc BEso1 nicely 
matches our benchmark TDDFT calculation at low and 
intermediate wavevectors, so we recommend using the 
model PBEsol xc hole, in all solid-state calculations 
where the hole is needed but full nonlocality is not impor- 
tant, instead of the more expensive TDDFT. We recall 
that the system-averaged hole, unlike the energy density, 
is uniquely defined, and is an observable at full coupling 
strength^. 

Figure [5] exhibits the wavevector-resolved correlation- 
only PBE, PBEsol, exact-RPA, and TDDFT surface en- 
ergies. We observe that at long wavelengths (k — > 0), 
where the LSDA is known to fail badly, and at short 
wavelengths (large k), where RPA is wrong, the GGA's 
under consideration are considerably close to our bench- 
mark TDDFT calculations. At intermediate wavevec- 
tors, however, GGA's cannot describe j c accurately, al- 
though PBEsol has been shown in Fig. [5] to give a very 
good description of j xc . This is due to a cancellation 
of the errors introduced at these wavevectors within the 
exchange and correlation contributions to j xc , which al- 
most cancel each other—. 
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IV. CONCLUSIONS 

We have constructed a PBEsol angle-averaged xc hole 
n xc (v, u) that satisfies known exact constraints and re- 
covers the recently reported PBEsol xc energy functional. 
Our construction of the PBEsol xc hole begins from and 
appropriately modifies the sharp cutoff correlation hole 
of Ref. [H and the smooth exchange hole of Ref. @. We 
also generalize [see Eq. (|7J|] the sharp cutoff procedure 
for the correlation hole to any GGA which has a positive 
gradient expansion coefficient. 

We have found that our PBEsol xc hole describes accu- 
rately the wavevector-resolved xc jellium surface energy 



for all values of the wavevector, thus providing support 
for the PBEsol GGA for solids and surfaces. 
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